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We investigate the low-energy electronic transport across grain boundaries in graphene ribbons 
and infinite flakes. Using the recursive Green’s function method, we calculate the electronic trans¬ 
mission across different types of grain boundaries in graphene ribbons. We show results for the 
charge density distribution and the current flow along the ribbon. We study linear defects at vari¬ 
ous angles with the ribbon direction, as well as overlaps of two monolayer ribbon domains forming 
a bilayer region. For a class of extended defect lines with periodicity 3, an analytic approach is de¬ 
veloped to study transport in infinite flakes. This class of extended grain boundaries is particularly 
interesting, since the K and K' Dirac points are superposed. 
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I. INTRODUCTION 

The synthesis of graphene by chemical vapor depo¬ 
sition (CVD) on metal surfaces is the most widely 
used method for producing graphene sheets.CVD 
graphene, as any other solid grown by CVD, is espe¬ 
cially prone to the formation of grain boundaries (GBs) 
and extended defect lines, which hinder its electronic 
properties. 

Graphene is being proposed for a variety of new elec¬ 
tronic devices.However, the required high-quality 
electrical properties are affected by the formation of poly¬ 
crystalline structures.These structures are practi¬ 
cally unavoidable by the growth methods known so far.^’^ 
As such, the scattering problem of an electron off a grain 
boundary (GB) becomes a theoretical and an experimen¬ 
tal relevant one.^^“^^ 

The sp^ bonding structure of carbon atoms in graphene 
gives rise to extended topological defects that are typ¬ 
ically composed of pentagonal, heptagonal, and octag¬ 
onal rings of carbon atoms, together with distorted 
hexagons.GBs are in general neither perfect straight 
lines nor periodic, intercepting each other at random an¬ 
gles. However, periodic straight GBs and defect lines 
can also be observed in graphene,^ and more interest¬ 
ingly they can be controllably synthesized at precise lo¬ 
cations and orientations,lifting the prospects for the 
engineering of arrays of such defects that would allow 
us to manipulate the electronic valley degree-freedom in 
graphene. 

GBs are known to strongly influence the properties 
of graphene, namely its chemical, mechanical, and elec¬ 
tronic ones.^^^^^ GBs are expected to present different 
degrees of transparency to electron transport, depending 
on their microscopic details and on the relative orienta¬ 
tions of the grains separated by In 

fact, measurements of electronic mobilities of different 
CVD samples, have shown that their electronic proper¬ 
ties strongly depend on the details of the CVD-growth 


recipes.Interestingly enough, and of direct rele¬ 
vance to our work, recent research has probed the electric 
properties of single GBs.^^’^^’^^ 

In a revealing work, Yazyev et al?^ have studied elec¬ 
tronic scattering from a wide variety of periodic GBs. In 
that work, based on momentum conservation along the 
periodic grain boundary, the authors have shed light on 
whether low-energy electrons travelling from one grain to 
the other may feel a transport gap at the GB. Their con¬ 
clusions were also quantitatively corroborated by first- 
principles quantum transport calculations (based on den¬ 
sity functional theory and the non-equilibrium Green’s 
function formalism). 

Our approach to the scattering problem due to GB’s 
follows two different routes. In the first route we use the 
recursive Green’s functions method to numerically cal¬ 
culate the transmission through defect lines in graphene 
ribbons. We map charge density over each sublattice site 
and also the current density through the defects. Fol¬ 
lowing the experiments,we consider grain bound¬ 
aries composed of extended linear defects of type 585 
(pentagons and octagons) and 5757 (pentagons and hep¬ 
tagons). Our results for the resistance across a linear 
defect compare well to recent experimental results.We 
also consider graphene ribbons with bilayer GBs: a spa¬ 
tial region where the grain boundary is composed by the 
superposition of two monolayer domains, as shown in 
Fig. 1. For this kind of overlapping bilayer boundary, 
previous results have already shown interesting conduc¬ 
tance oscillations.^^ Here we show that the transmission 
through these superpositions is reduced in comparison to 
the transmission in the single-crystal domains. We also 
present a spatial map of the current and the charge dis¬ 
tribution through these overlapped regions, which helps 
in the understanding of the transport properties of these 
systems. 

In the second route, we will concentrate on a particu¬ 
lar class of extended grain boundaries briefly addressed 
by Yazyev et al, namely, those grain boundaries with 
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periodicities such that both Dirac points (at each side 
of the grain boundary) are mapped into the F point of 
the projected Brillouin zone. In such cases, there will be 
intervalley scattering of massless Dirac fermions at the 
grain boundary. 

We have chosen to investigate zigzag aligned extended 
linear defect lines separating two grains with the same 
orientation (also referred to in the literature as degener¬ 
ate, i. e., zero misorientation angle, grain boundaries). 
Several such defect lines were proposed in the context of 
ab-initio works both on graphene and on boron nitride: 
the t7t5 defect line^^ and the 7557 defect line^^ (see Fig. 
8) are two such defect lines. As we will see ahead, this 
defect line allows for valley scattering to occur and can 
thus be regarded as a useful nanostructure for valleytron- 
ics circuits. 

For this second type of problems, and in the context 
of graphene’s single particle first neighbor tight-binding 
model, we will use the transfer matrix formalism^^ to an¬ 
alytically compute the transmittance of electrons across 
these grain boundaries. The boundary condition seen by 
the electrons at the grain boundary will be determined 
from its microscopic tight-binding model. In doing so, we 
will follow the methodology developed for the cases of the 
pentagon-only, zz(558), and zz(5757) defect lines.We 



FIG. 1. (color online) Schematic representation of the differ¬ 
ent linear defect structures considered as grain boundaries. 
On the top is the 5757 structure. In the middle we show 
the 585 linear defect. At bottom is the bilayer graphene of 
length L, formed by the overlap of two monolayer regions. 
The shadow areas represent the left and right semi-infinite 
contacts. 


will see that the obtained boundary condition explicitly 
introduces intervalley scattering. 


II. TIGHT-BINDING MODEL AND 
TRANSPORT FORMALISM 

We consider graphene’s tight-binding Hamiltonian 
H = —ti^j (c|cj + H.c.) 

<hj> 

where q (cj ) annihilates (creates) an electron at site i and 
(i,j) stands for pairs of nearest-neighbor atoms. We use 
the value of tij = 2.7 eV for the inplane nearest-neighbor 
hopping parameter and, when modeling the region of bi¬ 
layer graphene, we use =0.381 eV for the interlayer 
coupling.The extended linear defects are constructed 
by rearranging the positions of the atoms in the defect 
region. As shown in Fig. 1, this modifies the topology of 
the lattice and thus changes pristine graphene’s Hamil¬ 
tonian in Eq. (1). 

Coherent transport across grain boundaries in 
graphene is studied within the Landauer-Biittiker for¬ 
malism, which relates the conductance G{E) at a given 
energy E to the transmission function T{E) between the 
contacts as 


G{E) = GoT{E ), (2) 


2 -1 

with Go = 2^ ^ context of the previously 

referred first approach to scattering problems by GBs, 
the transmittance T{E) is evaluated by means of the re¬ 
cursive Green’s function approach using a two-terminal 
device configuration with contacts represented by the 
semi-infinite ideal graphene leads 


T = Tr 




(3) 


where Gs is the retarded Green’s function of the system, 
given by^^ 


Gs = [E'I-Hs-^l-^r]-^. (4) 

In these expressions is the Hamiltonian for the scat¬ 
tering region, = t^gL{R) stand for the self-energies 

coupling the scattering region to the leads, while E' is a 
shorthand for E' = E-\-ir]^ with 77 0. The self-energies 

and the broadening function l{r) = ~ 

are calculated from the electrode’s Green’s function gL{R) 
also obtained numerically using a recursive technique. 

Charge and current are intimately related through the 
continuity equation. The connection with the Green’s 
function arises from the quantum statistical average of 
the bond charge current operator, Jij = ^ — 

which is related to the lesser Green’s function 
In a steady state the bond charge current 
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including spin degeneracy is: 




dE 


t^^G<{E)-UjG<{E) 


( 5 ) 


where E^ = Ep ± eV/2^ while Iq stands for the natural 
unit of bond charge current density being given by Iq = 
2elh ll.biiE/eV. 

The lesser Green’s function in the absence of in¬ 
teractions can be solved exactly giving G^{E) = 

Gs{E)[rLfL + TnfnlGliE), where fuR) is the Fermi 
distribution of the left (right) contact and tji is the hop¬ 
ping parameter between sites j and i. The bond current 
Jij can be visualized as a bundle of flow lines bunched 
together along a link joining the two sites. 

Complementary to the current density, the charge den¬ 
sity at site j can also be expressed using the lesser Green’s 
function as: 


pEF-\-eV j2 

= ( 6 ) 

JEf-cVG 


It is noteworthy that at low bias and low temperature 
the charge density p, has the same distribution of the 
local density of states (LDOS). Given that we are inter¬ 
ested in how charge and current distributions are related, 
to keep explanations and figures as simple as possible, we 
will refer from now on to LDOS as charge distribution, 
with no loss of generality. 

In addition to the Landauer-Buttiker formalism [see 
Eq. (2)], it has been shown^^dQ Pow can we compute 
the low-energy limit of the conductance across this kind 
of defect lines. Interestingly, at low temperatures, the 
conductance across a defect line of size W [see Fig. 1] 
turned out to be linear in KpW and proportional to the 
transmittance [see Eq. (3)] close to Dirac point (E 0): 


G(E) = W 


QvQs 

Ait 


E 

hvp 


GoT{E) . 


( 7 ) 


The gate voltage Vg is nothing more than the spatial 
potential distribution created by the substrate’s charge 
distribution. We have estimated Vg for the GBs from the 
capacitor law 


_ qnd 

eA 


( 8 ) 


where n stands for the carrier density, d is the thick¬ 
ness and e is the dielectric constant of the substrate. In 
order to convert the experimentally measured gate volt¬ 
age into carrier density n we use the relation Vg = n/a^ 
where a = 2.5 x 10^^ m~‘^V~^ is a geometry-related fac¬ 
tor. From here onward, we will only consider the carrier 
density to be n = gsgvKp/dir^ where Kp i^ the momen¬ 
tum at the Fermi energy and g^ (gs) stands for the valley 
(spin) degeneracy. Finally, as in graphene Ep and Kp 
are proportional at low-energy, Ep = hvpKp^ then 

Y (p\ _ l 9s9v 

^9\-^) — a Atv yhvF^ ’ 


From Eq. (7) and Eq. (9) one can expect that the resis¬ 
tance of a periodic defect line, at low temperatures and 
in the linear regime, should behave as 


R=-«. 


1 

W9' 


( 10 ) 


This square root dependence on Vg should be clear from 
the experimental measurements of the resistivity across 
a grain boundary. 


III. MODIFIED CONDUCTANCE 
QUANTIZATION IN THE PRESENCE OF THE 
LINEAR DEFECT 


The study of scattering by extended defects is becom¬ 
ing increasingly more relevant, specially after the recent 
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FIG. 2. (color online)(a) Resistivity across the grain bound¬ 
aries structures 5757 and 585, compared to the experimental 
results of Tsen et al}‘^ . Inset: Transmission function. Here 
we used an armchair ribbon with gs = 2, g^ = 1, A" = 56 
atoms along the width and making further correspondence 
to kF = l/xm. (b) Spatial distribution of charge densities and 
current densities over the GB for different energies: E = 10~^ 
eV, F = 7 X 10“^ eV and F = 14 x 10“^ eV . The charge 
densities are schematically represented here for a narrower 
ribbon. The current densities are evaluated at different sites. 
The color of the arrow represents the magnitude of the electric 
current between any two neighboring sites, which are linearly 
normalized to the maximum value. 


(9) 
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work by Tsen et where the authors made electric 

measurements across a single grain boundary. They have 
found that the transport properties of these systems are 
strongly dependent on the GB’s microscopic details. 

Using Eq. (7) and Eq. (10) we calculate the resistance 
across two linear defects, the 5757 and the 585, and com¬ 
pare them with the experimental result from Tsen et al. 
(black circles) - Eig. 2(a). We verify that they agree to 
a good extent for Vg not too close to the Dirac point. 
The disagreement (at low carrier densities) between the 
experiment and our prediction is due to the effect of pud¬ 
dles, which are not taken into account in our calculations, 
and dominate the bulk region of the device. Note that 
the effect of puddles is increasingly relevant when one 
approaches energies very close to the Dirac point. In or¬ 
der to obtain the above results, we have considered an 
armchair ribbon with N = 8132 sites (^^ = 2, = 1 

and W ~ 1/im). The electronic properties at low-energy 
regimes, are obtained by a rescaling of the electronic 
properties of an armchair ribbon with N = 500 sites. 

The transport properties of these different linear de¬ 
fects are qualitatively similar. In the inset of Eig. 2(a) 
we show the electronic transmission and how it depends 
on the detailed geometry of the GB. In the continuum 
low-energy limit, both the 585 and the 5757 defects have 
a metallic behavior with a flat band crossing the Eermi 
level. 

The spatial distribution of charge density (for different 
energies) of a 585 linear defect is shown in Eig. 2(b). 
The density on each atomic site is represented as a disk. 
The different colors (red and blue) identify the sublattice, 
while the magnitude of the disk’s radius is proportional 
to the charge density at that site. We can see that the 
closer we are to the Dirac point {E = 10“^eV) the more 
localized the charge is in the region of the GB. We also 
plot on Eig. 2(b) the distribution of charge density for 
higher energies. In this case, a higher dispersion of the 
charge is apparent mainly before the 585 structure: the 
line defect acts as a potential wall. 

The corresponding current densities are also shown in 
Eig. 2(b), being evaluated at different sites using Eq. (5). 
The color of the arrow represents the magnitude of the 
electric current between any two neighboring sites, which 
are linearly normalized to the maximum value, accord¬ 
ing to the greyscale bar. Eor all the plotted energies, we 
observe that before and after the linear defect, the cur¬ 
rent flows in a specific horizontal pattern along armchair 
paths (streamlines) skipping some horizontal bonds, in 
accordance with recent ab-initio calculations of the cur¬ 
rent densities in pristine armchair graphene ribbons. 
Here, with the linear defect, it is interesting to observe 
how the current gives priority to some paths within the 
defect line, also in a periodical pattern. One can also 
see that the current density is smaller for the first repre¬ 
sented energy (i. e.,E = I0“^eV),as would be expected 
due to the localized nature of the charge density around 
the defect for this energy. 


IV. LINEAR DEFECT ORIENTATION 


The transmission across linear defects is known to 
significantly depend on their orientation angle. 

In Fig. 3(a)-(d), we show schematic representations of 
four different orientations of a 585 extended defect in a 
graphene ribbon: 0 = 0°,30°,60° and 90°. In Fig. 3(e) 
we plot the transmission probability as a function of the 
energy for these different angles of incidence at the de¬ 
fect. For all cases, we see that the particle-hole symmetry 
is broken due to the translation symmetry breaking in¬ 
troduced by the defect. 




FIG. 3. (color online) Schematic representation of a linear 
defect 585 at different orientation angles (a) 0°, (b) 30°, (c) 
60° and (d) 90°. (e) Transmission probability for each angle. 
The ribbon has A = 56 sites in the width. 


For ^ = 0°, the presence of the defect line located 
in the middle of the ribbon does not alter the metal¬ 
lic character observed in the transmission spectrum of a 
pristine armchair ribbon of the same width.Fig. 3(e) 
also shows that for 0 = 30° and 60° there are regions of 
vanishing transmission (opening of transport gap) close 
to the Dirac point. To further investigate the origins of 
these oscillations in the transmission, in Fig. 4 we map 
the charge and current density distributions for the se¬ 
lected energies indicated (by the arrows I, 2 and 3) in 
Fig. 3(e). 

The first energy (arrow I) corresponds to a transmis¬ 
sion resonance at an energy E = —0.21 eV, which is 
typical of resonant tunneling structures, where the con- 
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FIG. 4. (color online) Spatial distribution of charge densities 
and current densities over the linear defect 585 with a ori¬ 
entation angle of ^ = 30°, to corresponding to the energies 
indicated by the numbered arrows in Fig. 3(b): E = —0.21 
eV, F; = 1 X 10“^ eV and E = 0.21 eV. 


tinuum background is strongly suppressed at the dis¬ 
crete state localization. The associated charge and cur¬ 
rent density distribution are shown in Fig. 4(a). A high 
current density exactly following the defect line reveals 
its metallic character. Close to the contacts and far 
from the defect line, the current flow splits again into 
streamlines.^^ A similar behavior (high current density 
located along the linear defect) is found for resonances 
at angle of 60° (not shown here). We also found that the 
resonances of the first channel in the angles 30° and 60° 
are robust structures independent from W. The present 
results suggest an image of Fano-type resonances, a fi¬ 
nite coupling between the localized state associated to 
the linear defect and the delocalized continuum states 
associated with the armchair ribbon. 

In Fig. 4(b), we show the charge and current den¬ 
sity distribution for the second selected energy [arrow 


2 in Fig. 3(e)], corresponding to a vanishing transmis¬ 
sion close to the Dirac point. Interference in charge and 
current density on both sides of the linear defect can be 
seen, as well as on the linear defect. For this low energy 
we observe that the charge distribution is highly local¬ 
ized in part of the defect. The current seems not to flow, 
with its maximum local values circulating around the oc¬ 
tagons where the charge is concentrated and around oc¬ 
tagons symmetrically positioned with respect to the mid¬ 
dle of the linear defect. Note that the directions of the 
small arrows representing the local current flow around 
the octagons is different for Fig. 4(b) (local loopings) 
and Fig. 4(a) (net flow). The backscattering is evident. 
Therefore we conclude that at low-energies, i. e., close to 
Dirac point, the suppression of the transmission, at the 
first electron-like plateau, is due to charge localization 
and backscattering of a defect-related mode of the 585 
defect line. 

Fig. 4(c) corresponds to an energy value of E = 0.21 
eV. One can see that the current throughout the ribbon 
is not uniform and forms ambiguous paths. At this en¬ 
ergy, the current-density amplitude also flows across the 
defect, being greater in the edge region than in the center 
of the ribbon. The flow is mostly perpendicular to the 
linear defect. Note that an electron can travel between 
the source and the drain via many different transport 
channels. The local electric current profile at a given 
energy is nothing more than the result from the inter¬ 
ference between all the active transport channels at that 
energy. In particular, the existence of current loops for 
some defect orientations (see panels of Fig. 4), simply 
results from particular interference patterns arising from 
the different blockade of distinct transport channels by 
the linear defect. Moreover, the particular local current 
patterns that are observed result from the interplay be¬ 
tween the different components of the nanostructure: lin¬ 
ear defect topology and orientation, edges type and width 
of the graphene ribbon. In Fig. 3(b), for 0 = 90°, the 
first plateau does not present interference oscillations, the 
transmission is reduced in the vicinity of the Dirac point, 
due to the coupling of extended states at the edges with 
localized ones at the defect line. The oscillations at high 
energy range are simply Fabry-Perot interference effects. 


V. BILAYER GRAPHENE AS A GRAIN 
BOUNDARY 

In this section we investigate the electronic transport 
properties of a grain boundary defined by an overlap 
between two semi-infinite monolayer graphene regions, 
forming a bilayer region as represented in Fig. 5(a). Such 
overlaps have been experimentally observed.Here, to 
focus on the effects of the bilayer region on the transmis¬ 
sion, we consider periodical boundary conditions, avoid¬ 
ing edge localization effects. 

Figure 5(b) and (c) shows the transmission T(E) as a 
function of energy across a bilayer region of length L cor- 
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responding to 80 and 320 atoms superposed, respectively. 
We consider both the AB (Bernal) and AA stacking cases 
for the overlap regions. For both of them oscillations in 
the transmission are observed, with their frequency in¬ 
creasing for increasing overlap length L, in agreement 
with previous calculations for similar overlaps.This 
can be qualitatively understood by remembering that 
in such systems, the transmittance is set by the wave- 
function matching at the monolayer-bilayer interfaces. A 
monolayer eigenstate incoming from the left is going to 
be partially transmitted into the bilayer region and par¬ 
tially reflected back into the left monolayer. The portion 
of the wave-function transmitted into the bilayer region 
is going to propagate (acquiring a complex phase) un¬ 
til the second interface (bilayer-monolayer) and there it 
will undergo a similar scattering process: it is going to be 
partially reflected back to the bilayer and partially trans¬ 
mitted into the right monolayer. The resultant standing 
wave, in particular, the weight associated with each of its 
components (channels), is going to be the direct result of 
this interference process and will thus strongly depend on 
the length of the bilayer region and on the wave-number 
associated with each of those channels. The phases ac¬ 
quired by each of the wave-function’s components of the 
bilayer region at the second interface are going to be 
smaller for shorter bilayer regions (i. e., shorter L). In 
such cases, oscillations in the transmittance will require 
a greater change of the eigenstates’ wave-number, i. e., a 
greater increase in energy, as observed in both panels of 
Fig. 5. For comparison, we also show in Fig. 5(b) and (c) 
the transmission through a pristine monolayer and pris¬ 
tine bilayers AA and AB (considering the same width, 
40 atoms, to which periodic boundary conditions are ap¬ 
plied). One can see that the transmission through the 
pristine monolayer and bilayers is always higher than the 
transmission throughout the overlapped region. This is 
due to the presence of the interfaces, that act as scat¬ 
tering centers decreasing the system’s transmission. Fig¬ 
ure 5(d) shows the band structure for pristine monolayer 
and pristine bilayers AA and AB (of same width and 
also with periodical boundary conditions), which helps 
in understanding the origin of the plateaus in the trans¬ 
mission for each case. At low energies, in the case of AB 
stacking there is only one conducting channel, whereas 
in AA stacking there are always two conducting channels 
for each valley. This partially explains why in general 
the low-energy transmission for the AB bilayer structure 
is smaller than that for the AA bilayer structure. But 
in addition to this, we can easily check that the bound¬ 
ary condition at a monolayer-bilayer A A interface can be 
completely satisfied at low energies without the need for 
reflected components in the monolayer region. The same 
does not happen for the case of the monolayer-bilayer AB 
interface. Therefore, the upper bound for the transmis¬ 
sion at low energies is smaller in the bilayer AB case than 
in the bilayer AA case. Nevertheless, and by appropri¬ 
ately choosing the bilayer region length, L, we can still 
make bilayer AB case’s low-energy transmission higher 


than that of the bilayer AA case. 




FIG. 5. (color online) (a) Schematic representation of the 
overlap between two graphene monolayers, forming a bilayer 
region of length L. (b) Transmission throughout a bilayer re¬ 
gion of length L corresponding to 80 atoms superposed. Both 
AA and AB stackings are considered for the overlap (bilayer 
region). For comparison, it is shown the transmission through 
a pristine monolayer and pristine bilayers AA and AB of same 
width, (c) Same as in (b), showing now the transmission 
throughout a bilayer region of length L corresponding to 320 
atoms superposed, (d) Band structures for pristine monolayer 
and pristine bilayers AA and AB of same width. 


In Fig. 6 we map the spatial distribution of charge den¬ 
sity and local current density on each of the overlapping 
(AB stacking) graphene monolayer ribbons. Figure 6 (a) 
and (b) show such maps corresponding to the energies 
E = 0.001 eV and E = 0.1 eV, indicated by arrows 1 
and 2 in Fig. 5, respectively a minimum and a maxi- 
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mum values of trausmissiou iu the low-euergy limit. For 
E = 0.001 eV, i. e. iu Fig. 6(a), the charge deusity is lo¬ 
calized maiuly ou the left semi-iufiuite mouolayer, which 
correspouds to the bottom layer iu the bilayer regiou. Ou 
right semi-iufiuite mouolayer (i. e., the top layer of the 
bilayer regiou), the charge deusity rapidly decreases from 
the left to the right. Iu particular, its charge deusity is 
iusiguificaut iu the mouolayer regiou. Iu the regious with 
uou-vauishiug charge deusity of each layer, the charge 
deusity typically couceutrates ou ouly oue sublattice: the 
uou-dimer sublattice. This effect comes from the sublat¬ 
tice asymmetry iutroduced by the AB-stackiug iu bilayer 
grapheue.^^’^^ Moreover, for this euergy, the curreut deu¬ 
sity is very low ou both layers. Ou the other haud, for 
£■ = 0.1 eV (which correspouds to a peak iu the traus¬ 
missiou), we cau see iu Fig. 6(b) that agaiu the charge 
seems to be polarized ou the uou-dimer sublattice iu the 
ceutral part of the bilayer regiou. However, a careful ob- 
servatiou of other parts of the bilayer regiou shows that 
the charge is more homogeueously distributed over both 
sublattices there, thus allowiug electrou hoppiug betweeu 
sites aud betweeu layers,as observed iu the patteru 
of the zigzag curreut deusity streamliues. 

Similarly, iu Fig. 7 we show the charge aud curreut deu¬ 
sity distributiou for the case of the two mouolayers with 
au AA-stackiug overlap regiou. Figure 7(a) correspouds 
to E = 0.11 eV (the miuimum iu trausmissiou iudicated 
by the arrow 3 iu Fig. 5), while Fig. 7(b) correspouds 
to E = 0.2 eV (the resouauce iu trausmissiou iudicated 
by the arrow 4 iu Fig. 5). For both euergies, there is 
a clear charge wave aloug the bilayer leugth with charge 
oscillatiug betweeu the two layers. Similarly, curreut also 
oscillates betweeu layers. However compariug charge aud 
curreut deusities iu each layer, oue cau see au iuterestiug 
behavior: for the euergy correspoudiug to low trausmis¬ 
siou, Fig. 7(a) shows that there is a clear imbalauce, siuce 
the electrouic charge deusity aud the curreut deusity are 
couceutrated ou differeut parts of the bilayer regiou. Ou 
the other haud, iu Fig. 7(b), i. e., for the euergy corre¬ 
spoudiug to high trausmissiou, we agaiu observe a charge 
aud curreut oscillatiou betweeu layers aloug the leugth of 
the bilayer regiou, but uow these are iu phase, with the 
maximum curreut deusity spatially coiucidiug with the 
maximum charge deusity. Also uote the higher curreut 
deusities associated with the later euergy (see differeut 
curreut grey scale bar). 


VI. TRANSPORT ACROSS THREE PERIODIC 
GRAIN BOUNDARIES: TRANSFER MATRIX 
APPROACH 

Iu this sectiou we use the trausfer matrix formalism^^ 
to study the electrouic trausport across exteuded graiu 
bouudaries iu the scope of the siugle particle first ueigh- 
bor tight-biudiug model of (iufiuite) mouolayer grapheue. 
This method is eutirely equivaleut to the recursive 
Greeu’s fuuctiou uumerical method used up to this poiut 


of the mauuscript. It reduces the electrouic scatteriug 
problem to a set of matrix mauipulatious easy to work 
out by auy computatioual algebraic calculator, aud thus 
give rise to au aualytic solutiou of the problem. 

We will couceutrate ou a particular class of zigzag 
aligued periodic GBs that separate graius with the same 
orieutatiou (also kuowu as zero misorieutatiou augle 
GBs) aud have a periodicity that allows for iutervalley 
scatteriug of low-euergy electrous. Wheu the periodicity 
of the zigzag aligued GBs is a multiple of 3, both Dirac 
poiuts (as well as the F-poiut) are mapped iuto kxCi = 0 
- see pauels (c) aud (d) of Fig. 9. Therefore, aud iu 
coutrast with what happeus for the peutagou-ouly, 585 
aud 5757 liuear momeutum couserva- 

tiou does uot forbid low-euergy electrous from scatteriug 
betweeu valleys. Nevertheless, aud if we waut to kuow 
how much iutervalley scatteriug is a particular GB goiug 
to geuerate, we ueed to explicitly compute the bouud- 
ary couditiou matrix origiuatiug from its tight-biudiug 
microscopic model. 

From the diversity of GBs belougiug to this class, 
we have choseu to iuvestigate some that have beeu re- 
ceutly suggested iu the coutext of ab-iuitio works both 
ou grapheue aud ou borou uitride: the 7557 graiu 
bouudary^^ aud the t7t5 graiu bouudary^^ (see Fig. 8). 
As we will see ahead, these GBs coutrollably scatter elec¬ 
trous from oue valley to the other aud cau thus be re¬ 
garded as a useful uauostructure for valleytrouics cir¬ 
cuits, wheuever the desire is to destroy valley polariza- 
tiou. 

A periodic GB preserves the crystal’s trauslatiou iu- 
variauce aloug the GB directiou. Therefore, by Fourier 
trausformiug the system’s tight-biudiug Hamiltouiau 
aloug this directiou, we cau cast the problem of elec¬ 
trouic trausport iu a 2D crystal with a GB at its ceuter, 
as a problem of electrouic trausport ou a quasi-ID crystal 
with a localized defect at its ceuter. 

Iu order to work out this problem we will proceed as 
iu Ref. 18. From pristiue grapheue’s Hamiltouiau we 
start by writiug the tight-biudiug equatious away from 
the graiu bouudary (see Figs. 12 aud 13 for uotatiou 
clarificatiou) 

- ^A(n) = W"JiB(n-l) + B(n), (11a) 

-jB{n) = A{n) + WAA{n + l), (11b) 

where we have used the uotatiou Z(n) = 
[Zi{n), Z 2 {n), Zs{n)]^ for Z = A, B (sub-lattice 
ideutifier). Note that the uotatiou is hidiug the depeu- 
deucy of the Ai aud Bi ou the momeutum aloug 
the graiu bouudary directiou. Iu the above expressious, 
e aud t staud respectively for the euergy aud pristiue 
grapheue’s hoppiug parameter, while n gives the uumber 
of direct lattice vector U 2 = (—I,\/3)a/2 trauslatious 
away from the defect (see Figs. 12 aud 13). The matrix 
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FIG. 6. (color online) Charge and current density distributions for the nanostructure composed of two graphene (monolayer) 
ribbons that partially overlapped (AB-stacked in the overlap/bilayer region). The two panels stand for two different energies: 
(a) E = 0.001 eV, which corresponds to a minimum of transmission (indicated by the arrow 1 in Fig. 5). (b) = 0.1 eV, which 

corresponds to a maximum of transmission (indicated by the arrow 2 in Fig. 5). 


Wa reads 


Wa 


1 1 0 
0 1 1 

gSi/c^a Q I 


( 12 ) 


We can write the above tight-binding equations in the 
form 



■ B(n) ■ 

— (TD-, 

A(n) 


. Ein) _ 

— VI 

B(n- 1) 

A(n) 

B(n-l) 

= Q 2 

1 1 

1 1 


where the matrices Qi and 


read 




0 


fl3 

-I3 

KWa) 

-I3 


{Wa) 

0 


(13a) 

(13b) 


(14a) 

(14b) 


with I 3 standing for the 3x3 unit matrix. 

Eqs. (13) can be written in the form of a transfer 
matrix equation^^’^^ relating amplitudes at the atoms of 
the unit cell located at (n — l)u 2 with the amplitudes 
at the atoms of the unit cell located at nu 2 . Such an 
equation reads 


L(n) = T(e, kx)'L{n - 1), (15) 


withL(n) = [Ai{n),Bi{n)A 2 {n),B 2 {n)A‘i{A^Bz{n)Y', 
and the transfer matrix, T(e, given by 

T(e,A:^) = i?.Qi.Q2.i?^. (16) 

In the above equation, matrix R is simply used to 
change from the basis {^i, ^ 2 , ^ 3 , Ai, A 2 , A 3 } to the ba¬ 
sis {Ai, 5i, A 2 , ^ 2 , A 3 , ^ 3 }. It is written in Eq. (Al). 

Eollowing the method used for the cases of the 585 
and 5757 defect lines,we can find a basis where the 
transfer matrix becomes block diagonal with three 2 x 2 
matrices on its diagonal. In this basis the three modes 
of the problem are uncoupled. Moreover, around kx = 0 
two of these modes are low-energy (corresponding to each 
of the two Dirac cones), while the other is a high-energy 
mode. 

We can understand this fact from Eig. 9(b) where 
one represents the Eirst Brillouin zone (EBZ) originating 
from a honeycomb lattice whose direct vectors are chosen 
to be 3ui and U 2 . In such a EBZ, the two Dirac points 
are located at the same value of kx^ i. e., at kx = 0. It is 
thus natural that when setting kx = 0 in the transfer ma¬ 
trix given by Eq. (16), one obtains a transfer matrix that 
describes simultaneously low-energy electrons at each of 
the two valleys (together with an additional high-energy 
mode associated with the F-point region of the spectrum 
of pristine graphene). 
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FIG. 7. (color online) Charge and current density distributions for the nanostructure composed of two graphene (monolayer) 
ribbons that partially overlapped (AA-stacked in the overlap/bilayer region). The two panels stand for two different energies: 
(a) E = 0.11 eV, which corresponds to a low transmission (indicated by the arrow 3 in Fig. 5). (b) = 0.2 eV, which 

corresponds to a maximum of transmission (indicated by the arrow 4 in Fig. 5). 




FIG. 8. Scheme of two grain boundaries proposed in the con¬ 
text of ab-initio works both on graphene and on boron ni¬ 
tride. (a) the 7557 defect line^^ and (b) the t7t5 defect line.^^ 
In these schemes we highlight in blue the region of the grain 
boundary. 


0, while /+ and I— stand for the two modes with low- 
energy when kx ^ one associated with the K_ valley 
and the other with the K+ valley. The matrix A{kx) is 
explicitly written in Eq. (A2). 

As previously stated, in this basis the transfer matrix, 
T(e, kx)^ is block diagonal and reads 


T(e, kx) 


T?,(e, fca;) 0 0 

0 Ti-{e,k^) 0 

0 0 Ti+{e,k^) 


,( 19 ) 


where the three transfer matrices associated with each of 
the uncoupled modes are written in Eqs. (A3). 


The symbol A{kx) stands for the matrix mediating the 
change to the basis uncoupling the modes of the transfer 
matrix 

L(n) = A(A),)L(n). (17) 

We will denote the states in this new basis as 

L = [Ah, Bh, Ai_, Bi_, Aij^, (18) 

with h identifying the mode with high-energy when kx ~ 


A. The transmittance across the 7557 and the t7t5 
grain boundaries 

In a similar manner, the tight-binding Hamiltonian de¬ 
scribing the electronic structure close to the grain bound¬ 
ary can be used to write the tight-binding equations for 
the defect. With these we can compute a boundary con¬ 
dition relating amplitudes on either side of the defect 

L(l) =M.L(-1), 


( 20 ) 
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FIG. 9. (a) FEZ of pristine graphene (whose direct lat¬ 

tice vectors can be chosen to be ui = a(l,0) and U 2 = 
a(—1, \/3)/2). (b) FEZ of graphene whose unit cell has thrice 
the size of that of pristine graphene in the zigzag direction 
(direction of ui). In this case, the direct lattice vectors can 
be chosen to be 3ui and U 2 (as done in Figs. 12-13). (c) 
Spectrum of pristine graphene projected along the kx di¬ 
rection (parallel to the GE). (d) Gomparison between the 
-projected spectrum of (unfolded) pristine graphene (light 
green) and the same spectrum after a triple folding (red): in 
the latter, the two valleys (and the F-point) are mapped into 
kxa = 0 . 


where the boundary condition matrix, M, is a 6 x 6 ma¬ 
trix that depends both on the energy e, x-momentum 
kx and the electron hoppings characteristic of the grain 
boundary. In Appendix B we compute these matrices for 
the two grain boundaries we are investigating: the 7557 
and the t7t5 grain boundaries (see Figs. 12 and 13). 

Note that by expressing this boundary condition ma¬ 
trix in the basis that uncouples the modes of the transfer 
matrix T{e^kx) 


M = A(/ca^a).M.[A(/ca^a)] (21) 


we can conclude that in general, it mixes all the three 
modes of matrix T(e, kx). 

Given this, we now have all the ingredients needed to 
compute the coefficients involved in the electronic scat¬ 
tering by such defects. An incoming electronic wave from 
n = — oo, will be scattered by the defect at n = 0 produc¬ 
ing a reflected and a transmitted component. The wave- 
function on each side of the defect will then be given by 


L(n < 0) = , (22a) 

i=i 

L(n > 0) = ^ A”-i , (22b) 

i=i 

where pij and Tij are, respectively, the reflection and 
transmission scattering amplitudes from an incoming 
(from n = —oc) state, , into reflected, and trans¬ 
mitted, outgoing states. Finally, by imposing the 
corresponding boundary condition [see Eqs. (20) and 
( 21 )], we can compute the coefficients pij and Tij for a 
given energy and a given longitudinal momentum. 

For both the t7t5 and the 7557 grain boundary we 
have set the hopping parameters in the region of the 
grain boundary by estimating the corresponding carbon- 
carbon distances originating from the ab-initio results of 
Refs. 32 and 33, and then using the parametrization^^ 

exp[-a3 X (r“/ - (23) 

no 

where Vij stands for the distance between the carbons 
labeled by i and j (given in units of angstroms), the 
adimensional parameters 0^2 = 1.2785, 0^3 = 0.1383, = 

3.4490, while uq is the carbon-carbon distance in pristine 
graphene (in units of angstroms). 

In Fig. 10 we present the transmission probability for 
the 7557 grain boundary (see scheme of Fig. 12 and 
Appendix B 1 ) of an incoming electron of the K+ val¬ 
ley. The several transmittance curves of this flgure cor¬ 
respond to different energies and were drawn using the 
following hopping parameters at the defect: ^ = 0.98, 
7 = 0.94 and P = 0 . 1 . One can see on the several panels 
of this flgure that the intervalley scattering is comparable 
to the valley-preserving scattering. Both of them strongly 
depend on the energy and incidence angle, mainly due to 
the dependence on energy and kx of the boundary con¬ 
dition matrix (see its computation in Appendix B 1 ). 

Similar plots are presented in Fig. 11 for the t7t5 
grain boundary (see scheme of Fig. 13 and Appendix 
B2). These were obtained with the following hopping 
parameters: = 1.06, ^2 = 0.95, ^3 = 0.83, <^4 = 0.80, 

^5 = 1.30, ^6 = 1-05, 7 i = 1.23, 72 = 1.20, 73 = 1.18 and 
7 i = 1.36. In general, this choice of hopping parameters 
gives rise to a lower intervalley scattering at low energies 
than what is obtained for the 7557 grain boundary. 

We finalize by noting that the above transmittance 
curves are strongly dependent on the choice of the hop¬ 
ping parameters. In particular, at low energies these are 
the only parameters determining the boundary condi¬ 
tion matrix and therefore controlling the system’s trans¬ 
parency to incident electrons. Moreover, the GB’s scat¬ 
tering profile can be strongly enhanced or suppressed by 
small changes of the GB’s hopping parameters. There¬ 
fore, we may expect that intervalley scattering at the GB 
is deeply sensitive to modifications of the lattice’s geom¬ 
etry (namely strain) in the vicinity of the GB. 
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FIG. 10. Transmittance in terms of the incidence angle (in¬ 
coming electron chosen to be on the valley K+ Dirac point) 
for the 7557 grain boundary. The panels (a), (b), (c) and 
(d) correspond to scattering processes occurring at energies 
of, respectively, O.Olt, O.lt, 0.3t and 0.5t. The hopping pa¬ 
rameters (see Fig. 12) were set at: ^ = 0.98, 7 = 0.94 and 
/3 = 0.1. The dark blue curve stands for the transmittance 
preserving the valley degree of freedom (electron from the K+ 
below the GB scattering scatters to the same valley above the 
GB, i. e. T++ = |t++|^). The dashed light blue curve stands 
for the intervalley transmittance (electron at the K+ valley 
below the GB scattering into the valley K_ above the GB, i. 

e. T_+ = {v-/vj^) \t _ ^1^; v± stands for the velocity of the 

mode /±). 
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Appendix A: The bulk tight-binding equations 


The matrix R changing from the basis 
{Bi (n ), B 2 (n), Bs (n), Ai (n), A 2 (n), A 3 (n)} into the 
basis {Ai{n), Bi{n), A 2 {n)^ B 2 {n), As{n), Bs{n)} reads 


"OOOIOO' 

1 0 0 0 0 0 
0 0 0 0 1 0 
010000 • 

0 0 0 0 0 1 

0 0 1 0 0 0 

0/71 d/K 


(Al) 



FIG. 11. Transmittance in terms of the incidence angle (in¬ 
coming electron chosen to be on the valley K+ Dirac point) 
for the t7t5 grain boundary. The panels (a), (b), (c) and 
(d) correspond to scattering processes occurring at energies 
of, respectively, O.Olt, O.lt, 0.3t and 0.5t. The hopping pa¬ 
rameters (see Fig. 13) were set at: = 1.06, ^2 = 0.95, 

^3 = 0.83, ^4 = 0.80, ^5 = 1.30, ^6 = 1.05, 71 = 1.23, 
72 = 1.20, 73 = 1.18 and 71 = 1.36. The dark blue curve 
stands for the transmittance preserving the valley degree of 
freedom (electron from the K+ below the GB scattering scat¬ 
ters to the same valley above the GB, i. e. T++ = k++n. 
The dashed light blue curve stands for the intervalley trans¬ 
mittance (electron at the K+ valley below the GB scattering 
into the valley K_ above the GB ,i. e. r_+ = («_/«+) |r_+p; 
v± stands for the velocity of the mode /±). 


The matrix mediating the basis change that uncouples 
the modes of the transfer matrix, A{kx)^ reads 


A 



1 


1 0 
0 1 
1 0 
0 1 
1 0 
0 1 


g-70-27r/3)-^ 

l_|_gi7r/3 

0 

_g-i((/)-7r/3) 

0 

_g-d0+7r/3) 

0 


0 

g-i(0-2^/3)-^ 

l+e^^/3 

0 

_g-i((/)-7r/3) 

0 

_g-d0+7r/3) 


l^g-z7r/3 

0 

_g-i(20+7r/3) 

0 

_g-i(2(/)-7r/3) 

0 


0 

e-i 2 { 4 >- 2 ./ 3 ').^ 

l^g-i7r/3 

0 

_g-i(20+7r/3) 

0 

_g-i(2(/)-7r/3) 


(A 2 ) 


where = kxa. 


As just said, in this basis the transfer matrix, Eq. (16), 
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becomes block diagonal with three 2 x 2 matrices in its 
diagonal. The three pairs of modes, h, /+ and , de¬ 
couple and propagate independently. If we put ourselves 
around the Dirac point K+ = (0, — l)z^ 47 r/( 3 \/ 3 a), the 
upper matrix corresponds to the high-energy mode, the 
middle one corresponds to the Dirac cone identified by 
z/ = — 1 , while the lower matrix stands for the cone iden¬ 
tified by V = -1-1. For a general energy and momentum 
these three matrices read 


Tft(e,(/)) = 


1 


1 + 

T,_(e,</.) = /(0) 


— 1 —e 

e — 2 — 2 cos ( 


e + 


2 I e T-1 


/( 0 ) 


Ti+(e,(/)) = g{(^) 


— 1 —e 

2 I ^-1 


, (A3a) 
(A3b) 

(A3c) 


where we have again used (j) = kxa and have defined /( 0 ) 
and g{(j)) as 


piTv/S _ p-icf) 

f{<P) = --, 

1 — 2 cos (p 

(A4a) 

p-iTTl3 _ p-i4> 


^ • 

1 — 2 cos (p 

(A4b) 


Appendix B: The boundary condition of the 7557 
and t7t5 grain boundaries 


In this appendix we will briefly compute the boundary 
condition matrix associated with the two grain bound¬ 
aries investigated in Section VI. 


1. The boundary condition of the 7557 grain 
boundary 


Let us start by computing the boundary condition ma¬ 
trix relating the wave-function amplitudes on either side 
of the 7557 grain boundary (see Fig. 12 for a scheme of 
its crystalline structure). The tight-binding equations at 
the defect region read 


--^B{0) = A{1) + X^r> + BA{0), (Bla) 
-^D = A’(A(0) + B(0))+C?D, (Bib) 
-^A(0) = Am + WJiB(-l) + BB(0), (Blc) 


n = 1 



n — 0 
n = —1 












FIG. 12. Crystalline structure of the 7557 grain boundary.^^ 
The region of the defect line is highlighted in blue. 


where we use the notation Z(n) = [Zi(n), ^ 2(^)5 ^ 3 (^)]^ 
with Z = A^B^D. The matrices A, B and G read 


A = 


^00 
0 0 ^ 
0 0 0 


B = 


0 0 0 

0 13 0 

0 0 0 


G 


0 0 
Q Q 

0 0 0 


(B2a) 


(B2b) 


(B2c) 


where 0 = kxa. 

The above equations give the boundary condition be¬ 
tween either side of the grain boundary which reads 


B(l) 

A(l) 


B(-l) 

A(-l) 


(B3) 


where the are 6 x 6 matrices that read 


-V 


f2 

-Q 


^2 — 


- W\) Q-^VW^ 


-1 


-IWa 


-1 




(B4a) 

.(B4b) 


In Eqs. (B4) we have used the following definitions for 
the matrices V and Q, 


p = f I 3 - ^ (B5a) 

Q = B + X'^Tl , (B5b) 


where I 3 stands for the 3x3 identity matrix, while the 
matrix IZ reads 


7^ = - 


1 




0 -V 
0 


(B 6 ) 




































13 



1 

0 

-1 




A_ 




A 








FIG. 13. Crystalline structure of the t7t5 defect line.^^ The 
region of the defect line is highlighted in blue. 


Note that the above matrices depend on the reduced en¬ 
ergy, e/t, the longitudinal momentum, and the hop¬ 
ping parameters at the defect, 7 and /3. Similarly, 
the boundary condition connecting the two sides of the 
defect, i. e., L(l) = M 7557 .L(— 1 ), reads 

M7557 = (B 7 ) 

and in general depends on 7 , /3, e/t and k^- 


2. The boundary condition of the t7t5 grain 
boundary 

In Fig. 13 we can see the scheme of the crystalline 
structure of a t7t5 grain boundary. In what follows we 
will compute the boundary condition matrix relating the 
wave-function amplitudes on either of its sides. For such 
a grain boundary, the tight-binding equations in the de¬ 
fect region read 


- ^B(l) = G+A( 1 ) + WaA( 2 ), 

(B 8 a) 

-^A(l) = G+B(l) + AiA(l) + A 2 D, 

(B 8 b) 

-^D = a|A(1) + 5D + A3B(-1), 

(B 8 c) 

-^B(-l) = G_A(- 1 ) + AjD + A 4 B(- 

-l),(B 8 d) 

-^A(-l)=G_B(-l) + wiB(-2), 

(B 8 e) 


where, once more we use the notation Z(n) = 
[Zi(n), ^ 2 ( 71 ), Z 3 (n)]^, now for Z = A^B^D. The G± 
are diagonal matrices that can be written as G+ = 
diag[7i,72,73] and G_ = diag[ 72 , 73 , 7 i], while the Xi 


matrices (with i = 1,2,3,4) read 




X 2 


^3 


6 

■ 0 

■?2 


0 1 0 

10 0 , 
0 0 0 

0 

0 0 


0 C 4 


0 

0 

0 


r 6 0 01 
0 ^4 0 
000 ’ 

- 0 ^5^6- 




0 0 

0 0 0 

gi30 Q Q 


where again <f> = k^a. Finally, S reads 


5 = 


r 0 6 0 0-1 
6 0 0 0 
0 0 0 6 
.0 060. 


(B 9 ) 

(BIO) 


(Bll) 


(B12) 


(B 13 ) 


We can rewrite the above equations in a more compact 
form that allows us to write the equation relating the 
amplitudes at each side of the grain boundary (i. e., 

those at n = 2 with those at n = —2) in the following 
way 


B(2) 

A(2) 


N1.N2.N3.N4.N5.N6 


B(-2) 

A(-2) 


(B14) 


where the matrices Nj are now 6x6 reading 


and 


Ni 


N 2 


N 3 


fl 3 (WaA 

-I3 0 J ’ 

I{Wa)~" {Wa)~"g+ 

-I3 0 

g;Ai g+^X2P-^X3 

-I3 0 


N4 


Ns 


Ne 


Q-1F2 Q-^GZ^ 


-I3 

0 


1 1 

11 

1 

-I3 

0 

. 




-I3 

0 

5 


where we have used the following definitions 


Fi = 73+V+A2P-AI, 
F2 = h3+Xjp-^X3+X4, 

P = -lh-s, 

Q = XJP-^XI, 


(B15a) 

(B15b) 

(B15c) 

(B15d) 

(B15e) 

(B15f) 

(B16a) 

(B16b) 

(B16c) 

(B16d) 
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The above matrices depend on the reduced energy, e/t, 
the longitudinal momentum, and the hopping pa¬ 
rameters at the defect, and 7 j (with i = 1,..., 6 and 
j = 1,..., 4). 

It is now straightforward to write the boundary con¬ 
dition connecting the two sides of the defect L(l) = 
Mt 7 t 5 .L(—1), where the boundary condition matrix, M, 


is a 6 X 6 matrix given by 

Mt7t5 = i?.Ni.N2.N3.N4.N5.N6.i?^, (B17) 

where, for the sake of simplicity of notation, we have 
omitted the dependence of the matrices Mtnb and Ni on 
e/t, kx, Ci and 
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